Anomalous Nuclear Quantum Effects in Ice 
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One striking anomaly of water ice has been largely neglected and never explained. Replacing 
hydrogen ( 1 H) by deuterium ( 2 H) causes ice to expand, whereas the "normal" isotope effect is 
volume contraction with increased mass. Furthermore, the anomaly increases with temperature T, 
even though a normal isotope shift should decrease with T and vanish when T is high enough to 
use classical nuclear motions. In this study, we show that these effects are very well described by ab 
initio density functional theory. Our theoretical modeling explains these anomalies, and allows us 
to predict and to experimentally confirm a counter effect, namely that replacement of 16 O by 18 O 
causes a normal lattice contraction. 

PACS numbers: 31.15.A-,31.70.Ks, 65.40.De, 71.15.Pd 
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Numerous recent studies address the contribution 
of zero-point nuclear quantum effects to the structures of 
ice and water. The issue is delicate, because of the pe- 
culiar electrostatic-covalent nature [8j of the hydrogen 
bond (Hbond) in water. It is well known that hydro- 
gen bonded materials show an anti-correlation Q effect 
between the OH covalent bond and the OH-0 Hbond. 
As the OH-0 distance diminishes, the OH covalent bond 
weakens (its length increases and vibrational frequency 
diminishes [HIE!]) while the OH-0 Hbond does the op- 
posite. Closely related is the surprising expansion of the 
ice (Ih) crystal [l2|, 13 1 when H2O is replaced by D2O. 
A similar effect is observed in water, where the molec- 
ular volume of H2O is smaller than that of D20[14( at 
all temperatures. Zero point expansion of crystal lattices 
is a well understood phenomenon, and is almost always 
greatest for the lightest isotopes. For example, the vol- 
ume of 20 Ne expands by 0.6% with respect to 22 Ne at 
T=0 This "normal" isotope effect corresponds to 

a ~12% zero-point expansion of 20 Ne relative to a hypo- 
thetical "classical" or "frozen" lattice 3, 17 1. Since H2O 
and Ne have similar molecular masses, one might expect 
similar effects. However, the volume of H2O at T = is 
~0.1% smaller than that of D 2 0, It has rarely 



been mentioned in the literature that this is opposite to 
the usual behavior, and no explanation has been offered. 

In this paper, we explain this effect as an interesting 
coupling between quantum nuclear motion and hydrogen 
bonding, that may be relevant also to the structure of liq- 
uid water. Our analysis shows that, despite the anoma- 
lous isotope effect, quantum ice actually has a volume 1% 
larger than it would have with classical nuclei. The effects 
are smaller than in Ne mostly because of delicate cancel- 
lations. We exploit these cancellations to make critical 
comparisons of: (i) quasiharmonic theory versus fully an- 
harmonic path-integral molecular dynamics (PIMD); (ii) 
ab initio forces versus flexible and polarizable empirical 



force fields (EFF); and (iii) various flavors of ab initio 
density-functional theory (DFT) exchange and correla- 
tion (XC) density functionals (DF) with and without in- 
clusion of van der Waals (vdW) interactions. We find: 

(i) quasiharmonic theory is satisfactory for this problem; 

(ii) present state of the art EFFs are not good enough 
to describe nuclear quantum effects in water; and (iii) 
all the DFs considered describe qualitatively the anoma- 
lous effects, although some versions perform better than 
others. 

Within the volume-dependent quasiharmonic approx- 
imation (QHA), the equilibrium volume V(T) is ob- 
tained by minimizing at each T the Helmholtz free energy 
F(V,T) [3 [3: 



F(V,T) = E (V) + 

\hhJ k (V) 
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where Eq(V) is the energy for classical (T = or frozen) 
nuclei, at the relaxed atomic coordinates for each vol- 
ume, ujk are the phonon frequencies, with k combining 
the branch index and the phonon wave vector within the 
Brillouin zone. Their volume dependence is linearized as: 
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where 7^ = — (Vq /uik){duik/ dV)v is the Griineisen pa- 
rameter of the mode, and Vq is the equilibrium volume 
of Eq(V). Wfe(Vo) and 7fc(Vo) are obtained by diagonaliz- 
ing the dynamical matrix, computed by finite differences 
from the atomic forces in a (3 x 3 x 3) supercell, at two 
volumes slightly below and above Vq. As shown in the 
supplementary information [20j | (SI), this linearization is 
an excellent approximation to the full QHA, in which the 
frequencies are calculated at each volume. It can be easily 
shown 20] that at T=0 the quasi-harmonic volume shift 
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with respect to the frozen lattice is isotope dependent, 
and proportional to the average (jk^k)- For increas- 
ing temperature, quantum effects become less relevant, 
and the shift converges to the classical result, where the 
change in volume is oc T and isotope independent. There- 
fore, for high enough T, the differences between isotopes 
must vanish. 
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FIG. 1. (Color online). Volume per molecule for different 
isotopes calculated using the q-TIP4P/F force field both with 
PIMD @| simulations and the QHA. 



Anharmonic effects, beyond the QHA, could change 
the V{T) curves calculated with this method. However, 
experimental Raman spectra of ice Ih, as a function of 
pressure and temperature Io|, 21, 22 1 show a relatively 
temperature- insensitive stretching band, with dw/dT < 
0.6 cm K . Such a small anharmonicity should not 
modify the QHA results. 

To test the validity of the QHA in ice Ih, Fig. [T] com- 
pares the QHA result for V(T) with our fully anharmonic 
PIMD simulations @|, using the q-TIP4P/F force field @ 
in both cases. Differences between the QHA and PIMD 
increase with temperature, but their overall agreement 
is satisfactory. In particular, with this force field both 
calculations predict a normal isotope effect @, which 
includes a lattice contraction when H is replaced by D 
and a convergence of H and D volumes with increasing 
T, contrary to the experimental result. To complement 
these results, we have repeated the calculations using the 
polarizable TTM3-F[U [H potential. This force field 
has recently been shown to outperform q-TIP4P /F when 
compared to experiments in PIMD simulations of H2 18 
and D 2 18 0. Results provided in TableUand in the SI [2(ij . 
show that this polarizable force field also fails to repro- 
duce the anomalous isotope effect. However, as it will 
be discussed below, TTM3-F improves over q-TIP4P/F, 
displaying a stronger anti-correlation effect. 

Our DFT simulations are performed using the siesta 
method 
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and improved version of the basis set used in Ref. [27 . 
Given the precision needed to obtain accurate Griineisen 
parameters, which involve numerical second derivatives 
of the forces, we have used highly converged grids for real- 
space and k-sampling integrations. All our residual forces 
were smaller than 10~ 3 eV/A. Details of the parameters 
and convergence tests are provided in SI 20] . 

Hexagonal ice Ih is the most common form of ice, where 
oxygens occupy the hexagonal wurtzite lattice sites. The 
two covalently-bondcd protons have six possible orienta- 
tions, but are constrained by Bernal-Fowler "ice rules" 
to have one proton per tetrahedral 0-0 bond. This 
structure is characterized by proton disorder [28(. Ice 



26j, with norm-conserving pseudopotentials 
and numerical atomic basis sets for the valence electrons. 
We use a basis of triple-C polarized orbitals, a longer 



XI is the proton-ordered phase of ice that is formed by 
a transition from hexagonal ice Ih at 72 K if catalyzed 
by KOH~ 29]. The lattice parameters are slightly dif- 
ferent from those of ice Ih, with the hexagonal symmetry 
weakly broken $GJ- We have performed calculations for 
three crystal structures: (i) The Bernal-Fowler structure, 
with 12 molecules per unit cell, a net dipole moment 
along the c axis, and proton order consistent with the 
hexagonal symmetry, (ii) a version of the ordered ice XI 
structure with 4 molecules per unit cell, constrained to 
be hexagonal, also with a net dipole moment along the c 
axis; and (iii) a 96-molecule proton-disordered structure 
with no net dipole moment [6j. T he oxygen lattice disor- 
der observed in experiments [31| , is of the same order of 
magnitude as what we observe in our PIMD simulations 
of the proton disordered structure 0, [32j . 

As the results obtained for cells (i) and (ii) are very 
close, in the following we will present results for struc- 
tures (ii) and (iii) . Table Q] shows the T = volumes 
for several isotope combinations of H, D, 16 0, and 18 0, 
for all the DFT functionals [331436T ] and force fields used 
in this study. Also included in this table is a 32 beads 
PIMD result for a single unit cell (r sampling) using the 
PBE-DF, compared to a calculation with the QHA for an 
identical system. The simulation was done at T = 200 
K. 

As the experimental study from Rottgcr et at did not 
consider H2 ls O, we have performed high resolution X- 
ray diffraction experiments of the three isotopes H2O, 
D2O and H 2 ls O. Results are presented in the table. Ex- 
periments were performed with samples of water frozen 
in thin-walled glass capillaries, open at both ends to ac- 
commodate the volume expansion upon freezing, or in 
flexible polymide tubes. In all cases Si powder within the 
sample promoted nucleation of multiple crystals and pro- 
vided an internal standard for precise determination of 
X-ray wavelength. X-ray measurements were performed 
at a temperature of 100 and 220K K at the X16C high- 
resolution powder diffractomctcr at the National Syn- 
chrotron Light Source, using X-rays of nominal wave- 
length 0.7 A. Diffraction peak widths were on the order 
of 0.02° FWHM. Independent samples were prepared by 
slowly (~ 0.05 K/sec) cooling and by quenching with liq- 
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uid Nj. Lattice parameters were determined from mea- 
sured positions of 21 ice diffraction peaks. Reported er- 
ror bars encompass selected measurements at different 
temperatures. A detailed experimental study including 
temperature dependence will be published in a future 
contribution. 

Except for vdW-DF rovPBE , all the other XC function- 
al predict an anomalous isotope effect at T = when H 
is replaced by D, in agreement with experiments. How- 
ever, the isotope effect has the normal sign for the O 
atom. Our experiments confirm this result, with a 0.06% 
volume contraction when 18 O replaces 16 O (T=100K). 
The comparison of structures (ii) and (iii) shows that 
the results are largely independent of the ordering of the 
protons. Agreement with experiments improves as Bril- 
louin zone integration is improved. With respect to their 
generalized gradient approximation (GGA) counterparts, 
vdW-DFs soften the afore-mentioned anticorrelation ef- 
fect, reducing the magnitude of the H— >D isotope shift, 
but have little effect on the O shift. When phonons from 
the full Brillouin zone are accounted, the vdW-DF rcvPBE 
fails to predict the isotope shift at T = 0. However 
the anomalous shift for this functional is recovered at 
T > 100 K 0. Overall, vdW-DF PBE has proven to 
be very robust for a variety of structural and dynami- 
cal properties of water (s, 35, 37|. It also gives our best 
lattice constant for ice at T—0. Therefore, in the fol- 
lowing we use this functional and the QHA to explore 
the volume of ice Ih in structure (ii) , including full Bril- 
louin zone phonon integration, as a function of isotope 
masses and temperature. However, the results, and in 
particular the anomalous isotope effect, are very robust 
and largely independent of the XC functional, harmonic 
approximation, system size, and proton ordering. 
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FIG. 2. (Color Online). Volume change V(T)/Vn 2 o(0) - 1, 
relative to H2O at T — 0, for different isotopes calculated 
using the QHA with the vdW-DF PBE functional. Also shown 
are the experimental results from Ref. Q3- 

Fig. [5] shows V(T) of ice Ih for standard isotope sub- 
stitutions of H and O. Experimentally, the anomalous 



H— s-D isotope effect increases from 0.09% at T — 10 K 
to 0.25% at T = 250 K, and this increase is reproduced 
by our calculations (from 0.16% to 0.32%). The classical 
volume becomes larger than any of the quantum results 
above 100 K. This implies that, for larger temperatures, a 
classical isobaric ab initio molecular dynamics simulation 
of ice will overestimate its volume, relative to a quantum 
PIMD simulation. 
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FIG. 3. (Color online). Top: Density of vibrational states for 
H2O, projected onto H and O atoms, for the ordered ice Ih 
structure, as obtained with vdW-DF PBE functional. Bottom: 
average Griineisen constants of the different modes. 

To analyze the complex isotope-dependent free energy 
surface, we plot in figure [3] the vibrational density of 
states, projected onto the two atomic species, as well as 
the average Griineisen constants of the different phonon 
branches. The volume change upon isotope substitution 
at T = is approximately proportional [lj, to the av- 
erage value (uiklk)- Negative 7fc's imply a softening of the 
modes with increasing pressure, favoring smaller volumes 
for lighter isotopes. Thus, the anomalous isotope effect, 
F(H 2 0) < F(D 2 0), is due to the negative 7^ of the cova- 
lent OH stretching modes (oj ~ 3100-3500 cm -1 ) which 
have more than 95% H weight. Although the likewise H- 
dominated librational modes (oj ~ 600-1000 cm -1 ) have 
a positive 7^, they are not enough to balance the vol- 
ume shrinking contribution of the stretching frequencies. 
In the case of the TTM3-F force held, the two contribu- 
tions are very small and almost completely cancel out, 
effectively producing a very harmonic ice crystal, with 
a tiny anomalous effect at high T. On the other hand, 
q-TIP4P/F largely underestimates the OH-O/OH anti- 
correlation, predicting a normal isotope effect for all tem- 
perature ranges. Long range torsional modes, represent- 
ing rigid tetrahedral rotations with oj « 50 cm -1 , also 
have 7a_, < 0. These modes are responsible for the nega- 
tive thermal expansion below 70 K. The bending modes, 
with oik = 1700 cm -1 and 7^, sa have little effect on 
volume. 
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TABLE I. DFT volume (in A 3 /molecule) for proton ordered (H-ordered) and proton disordered (H-disordered) ice Ih for 
different isotopes, obtained with the quasiharmonic approximation (QHA) or path integral (PIMD) simulations, k-mesh is 
the effective number of k points for sampling the 4-molecule hexagonal Brillouin zone in the phonon calculation (one for re- 
sampling). IS(A-B)= y^gj — 1, is the relative isotope shift for the exchange of isotope A by B. The exchange and correlation 

(XC) functionals are: PBE vdW-DF PBE H3, HH|, revPBE [H, and vdW-DF rcvPBE The force fields (FF)are q- 

T1P4P/F0] and TTM3-F[2^] V c i a is the volume for classical nuclei. Also shown are the experimental results from ref [l^] and 
the ones obtained in this work. Note they are at different temperatures. 
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Substitution of 16 by ls O affects mainly the low fre- 
quency modes, dominated by positive values of jk, pro- 
ducing a normal isotope effect. The temperature depen- 
dence of the isotopic O substitution is also normal, and 
the volume shift is 50% smaller at T=220 K than at 
T=l 00 K. Somewhat surprisingly, the net effect at T = 0, 
relative to classical nuclei, is dominated by quantum oxy- 
gen, resulting in a quantum volume ~1% larger, This 
small expansion (f times smaller than that of Ne) is a 
consequence of two competing anharmonicities: the con- 
traction effect of H-dominated stretching modes and the 
expansion effect of librational and translational modes. 
As T increases, the contribution of the stretching modes 
becomes dominant, causing the net quantum effect to 
change sign and to become anomalous above ~70 K. 
This dominance increases with T, making the volume 
shift four times larger at the melting temperature than 
at T=0. These results are not inconsistent with the re- 
quirement that, at high enough T the isotope shift is 
isotope independent, but we find that the convergence 
towards the classical limit starts at T>^900 K. 

Our results may also have significant implications for 
the understanding of nuclear quantum effects in liquid 
water (in which the anomalous isotope shift is experi- 
mentally larger than in ice 14{). PIMD simulations, us- 
ing the q-TIP4P/F and TTM3-F EFFs, produce a less 
structured liquid than classical MD simulations @, OJ ■ 
However, as we have seen, these EFFs do not repro- 
duce the anomalous isotope effect in ice, because they 



fail to describe accurately the derivatives of the frequen- 
cies, which govern the anharmonicities and the nuclear 
quantum effects in the structure and dynamics. This sug- 
gests that these models may be inadequate to reproduce 
some quantum effects in the liquid as well as in the solid. 
Therefore, the observed loss of structure in the liquid, for 
quantum vs classical nuclei, should be reanalyzed with an 
EFF that reproduces the anomalous quantum effects in 



In conclusion, we have shown that the anomalous nu- 
clear quantum effects on the volume of ice can be fully un- 
derstood using the quasi-harmonic approximation with 
density functional theory. The study fully explains a 
rare, seldom mentioned, property of ice which should be 
included in the list of water anomalies |38j |. as an exam- 
ple in which quantum effects are anomalous and increase 
with temperature. 
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